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Abstract 

The Rayleigh conjecture about convergence up to the boundary of 
the series representing the scattered field in the exterior of an obstacle 
D is widely used by engineers in applications. However this conjecture 
is false for some obstacles. AGR introduced the Modified Rayleigh Con- 
jecture (MRC), which is an exact mathematical result. In this paper we 
review the theoretical basis for the MRC method for 2D and 3D obstacle 
scattering problems, for static problems, and for scattering by periodic 
structures. We also present successful numerical algorithms based on the 
MRC for various scattering problems. The MRC method is easy to imple- 
ment for both simple and complex geometries. It is shown to be a viable 
alternative for other obstacle scattering methods. Various direct and in- 
verse scattering problems require finding global minima of functions of 
several variables. The Stability Index Method (SIM) combines stochastic 
and deterministic method to accomplish such a minimization. 

Key words: obstacle scattering, Modified Rayleigh Conjecture, Stability 
Index Method. 
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1 Introduction 

In this paper we review our recent results on the Modified Rayleigh Conjecture 
(MRC) method. The method is applied to multidimensional obstacle scattering 
problems, as well as to scattering by periodic structures. Also we discuss an 
application of the MRC to static problems, and preliminary results on inverse 
obstacle scattering by MRC. Numerical results illustrate the performance of var- 
ious MRC algorithms. The paper concludes with a presentation of the Stability 
Index Method (SIM) for global minimization. 
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The basic theoretical foundation of the method was developed in [27]. The 
MRC has the appeal of an easy implementation for obstacles of complicated 
geometry, e.g. having edges and corners. In our numerical experiments the 
method has shown itself to be a competitive alternative to the BIEM (boundary 
integral equations method), see [13]. Also, unlike the BIEM, one can apply the 
algorithm to different obstacles with very little additional effort. A similar 
method is discussed in [10] 

We formulate the obstacle scattering problem in a 3D setting with the Dirich- 
let boundary condition, but the method can also be used for the Neumann 
boundary condition, corresponding to acoustically hard obstacles, and the Robin 
boundary condition. 

Consider a bounded domain Del 3 , with a Lipschitz boundary S. Denote 
the exterior domain by D' = M. 3 \D. Let a, a' <G S 2 be unit vectors, where S 2 is 
the unit sphere in M 3 . 

The acoustic wave scattering problem by an acoustically soft obstacle D 
consists in finding the (unique) solution to the problem (1.1)-(1.2): 

(V 2 + k 2 ) u = in D', u = on S, (1.1) 



e ikr /]\ x 

u = u + A(a', a) h o - ) , r := \x\ — > oo, a := — . (1.2) 

r \rj r 

Here u :— e lka ' x is the incident field, v := u — u is the scattered field, A(a' , a) 
is called the scattering amplitude, its fc-dependence is not shown, k > is the 
wavenumber. The scattered field v is an outgoing solution of the Hclmholtz 
differential equation (1.1), that is, a solution which satisfies the radiation con- 
dition 

dv 
d\x\ 



lim 

Denote 



|x|=r 



ds = 0. (1.3) 



A e (a) := / A(a',a)Y e {a')da', (1.4) 
Js 2 

where Yi(a) are the orthonormal spherical harmonics, Yi = Y( m , — £ < m < I. 
Let a ball Br := {x : \x\ < R} contain the obstacle D. Let he(r) be the 

spherical Hankel functions, normalized so that hi{r) <~ ^— as r — > +oo. In the 
region r > R the solution to (1.1)-(1.2) is: 



u(x,a) =e ika - x +J2Ma)*e(x), */(ar) := Y e (a')h e (kr), r > R, a' = -, 

(1.5) 

where r = |x|, the sum includes the summation with respect to m, — £ < m < £, 
and At(a) are defined in (1.4), see [26]. 

The Rayleigh conjecture (RC) is: the series (1.5) converges up to the 
boundary S (originally RC dealt with periodic structures, gratings). This con- 
jecture is false for many obstacles, but is true for some ([3, 22, 28]). For example, 
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if n = 2 and D is an ellipse, then the series analogous to (1.5) converges in the 
region \x\ > a, where 2a is the distance between the foci of the ellipse [3]. In 
the engineering literature there are numerical algorithms based on the Raylcigh 
conjecture. These algorithms use projection methods and are reported to be 
unstable. Moreover, no error estimate has been obtained for such algorithms. 
These algorithms cannot converge for arbitrary obstacles, because the Rayleigh 
conjecture is false for some obstacles. 

Our aim is to give a formulation of a Modified Rayleigh Conjecture (MRC) 
which holds for any Lipschitz obstacle and can be used in numerical solution of 
direct and inverse scattering problems. In other words, while the MRC still has 
the word "conjecture" in its name, it is a proven mathematical result for the 
scattered field in the exterior domain D' . In contrast to algorithms based on the 
invalid Rayleigh Conjecture, the MRC-based algorithms, like the ones described 
here, converge and an error estimate for the approximate solution they yield has 
been obtained in [27] (see also [34], Chapter 12). This error estimate is sharp 
in the order e. 



2 Modified Rayleigh conjecture 

What we call the Modified Rayleigh Conjecture (MRC) is actually the following 
Theorems 2.1, see [27], and 2.3, see [16]. We denote by H™ C {D') the set of func- 
tions from the Sobolev space H m (D) for any compact strictly inner subdomain 
D of D', so that the distance from D to S is positive, dist(D, S) > 0. 

Theorem 2.1. Let v = u — uq be the scattered field, where u is the solution to 
(1.1)-(1.2). Then there exists a positive integer L = L(e) and the coefficients 
C£ = Ci(e), < I < L(e) such that 



(%). 

where 

(n). 

and 
where 



\\u + v e \\ L 2 {s) <e, (2.1) 
v t (x) = J2c e (e)* e (x). (2.2) 

£=0 



v\\l'(S) < e (2-3) 
= O(e),e->0, (2.4) 



III • III - II ' \\h™ c (D') + II ' \\l2(D>:{1+\x\)->) 

7 > 1 , m > is an arbitrary integer. 
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(in). 

a(e) -> A t , as 0, W, 
where Ai := At(a) is defined in (1.4). 

Proof. First, we prove item (i). Then we establish Lemma 2.2, and continue 
with the proof of (ii) and (iii) . 

(i) Without loss of generality we can assume that the origin is an interior 
point of the domain D. To establish (2.1) it is sufficient to show that 

H :=spm{i$ e (s) : < I < oo, s g S} = L 2 (S) . (2.5) 

Suppose that there exists p G L 2 (S), p ^ 0, such that p _L H in L 2 (S). Define 
the single-layer potential by 

I s -2/1 

where (is is the surface area element. Let U C D be a ball centered in the origin. 
Then the addition theorem for the fundamental solution implies that W{y) = 
for any y G U. 

By the unique continuation principle W = in D. In particular W = on 
the boundary S. Since W is an outgoing solution of (V 2 + k 2 )W = in D' with 
W = on 5, one concludes from the uniqueness of solutions to the Dirichlet 
problem in D 1 that W = in M 3 . Finally, the jump properties of the normal 
derivative of the single-layer potential imply that p = in L 2 (S). We have 
followed the argument from [33], p. 160. □ 

Lemma 2.2. Given g g L 2 (S), let w be the outgoing solution of the exterior 
Dirichlet problem (V 2 + k 2 )w = 0, inD' w ~ g on S. Then there exists a 
constant C > 0, independent of w, such that 

\\H\\<C\\g\\ L2(S ), (2-7) 

where ||| • ||| := || • (£)') + II • llL 2 (_D' : (i+|a;|)-T); 7 > 1, m > is an arbitrary 
integer, and H m is the Sobolev space. 

Proof. Let G be the Dirichlet Green's function of the Laplacian in D': 

(V 2 + k 2 ) G = -5{x - y) in D', G = on S, (2.8) 



lim 

r^oo 



o\x\ 



ds = 0. (2.9) 

l\x\=r "F| 

Let be the unit normal to 5 pointing into D' . By Green's formula one has 
w(x)= j^g(s) — (x,s)ds, xeD>. (2.10) 
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The estimate for the HJ^ c (D')-norm part of (2.7) follows from this representation 
and from the Cauchy inequality: 



\D^w(x)\ < \\g\\ L 2 {S ) 



dD ( / } G 



dN 



x,s) 



< c(x)\\g\\ L 2 {s) , 



where c(x) < c(d) for all x £ D' such that the distance dist(x, S) > d > 0. 

For the L 2 -weighted norm part of (2.7) let R > be such that D C B R = 
{x e R 3 : \x\ < R}. Let D' R = Br \ D, and S R be the boundary of B R . The 
estimate 



dG 



< 



l + \x\ 



\x\ > R, 



and formula (2.10) imply 



\\w\\ L 2 (Sr) < c|| ff || L 2 (5) , 



(2.11) 



(2.12) 



where here and in the sequel c and C denote various constants. Also, using the 
Cauchy inequality, formula (2.10), inequality (2.11) and the assumption 7 > 1, 
one gets 



HU 2 (M>-R;(l + M)-^) < chh^s) 



1 



(i + M 



,7+1 



L 2 (M>,R) 



<c\\g\\ L2{s) . (2.13) 



To get the estimate for ||w||l 2 (d' ) choose R such that k 2 is not a Dirichlet 
eigenvalue of -A in D' R . Then ([20]* p.189): 

II^IIh-(d^) < c[||(A+/c 2 )w|| wm -2 (D / 7) + ||w;|| ffm -o.5 (SR) + ||w|| H ™-o.5 (s) ]. (2.14) 

The space H in the first term of the right-hand side in (2.14) is different 
from the usual Sobolev space, but this term is equal to zero anyway because 
(A + k 2 )w = 0. 

Let to = 0.5 in (2.14). Then 



\\w\\h° -'(d> r ) < 4\\w\\ L 2 {Sa) + \\w\\ L 2 {s) }. 
Since w = g on 5, then (2.12) and (2.15) imply 
\\w\\l2(d> r ) < c\\g\\ L 2 {s) . 



(2.15) 

(2.16) 
□ 



Proof of Theorem 2.1, continued. 

(ii) Inequality (2.3) is the same as (2.1), since v — — u n on S. Estimate (2.4) 
follows from (2.3) and Lemma 2.2. 

(hi) Inequality (2.3) yields the convergence of v e to v in the norm |j • ||i,2(s). 
By (2.12) ||u e - v\\ L 2 (Sr) -» 0, as e -» 0. On S R one has v = EZoM^e 
and v e — J2e=d c e^i- Multiply v e (R,a') — v(R,a') by Yg(a'), integrate over 
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S 2 and then let e — ► 0. The result is (iii), and the proof of Theorem 2.1 is 
completed. □ 
The difference between RC and MRC is: (2.1) does not hold if one replaces v e 
by X^=o Aela^e, and lets L — > oo (instead of letting e — > 0). Indeed, the series 
X^Lo ^e( a )^ 't diverges at some points of the boundary for many obstacles. 
Note also that the coefficients in (2.2) depend on e, so (2.2) is not a partial sum 
of a series. 

For the Neumann boundary condition one minimizes 



d[u + Y,t=Q c ^t\ 



ON 



l 2 (S) 



with respect to eg, and obtains essentially the same results. 

According to Theorem 2.1 the computation of the outgoing solution to (1.1)- 
(1.2) is reduced to the approximation of the boundary values in (1.1) by the 
linear combinations of the functions \P g restricted to the boundary S. A di- 
rect implementation of the above algorithm is efficient for domains D not very 
different from a circle, e.g. for an ellipse with a small eccentricity but it fails 
for more complicated regions. The numerical difficulties happen because the 
spherical Hankel functions hi with large values of I are bigger than hi with 
small values of I by many orders of magnitude. A finite precision of numeri- 
cal computations makes it necessary to keep the values of L not too high, e.g. 
L < 20. This restriction can be remedied by the following modification of the 
above algorithm, see [13, 16]: 

Theorem 2.3. Let v := u — u , where u is the solution to (1.1)-(1.2). Let 
e > 0, and L be a nonnegative integer. Suppose U is an open subset of D. 

Then there exist a finite subset {z\, z%, zj} C U, and the coefficients 
ce(e,Zj), < i < L, 1 < j < J = J(e), such that the following inequalities 
(2.17) and (2.20) hold: 

(i) . 

IK + v e \\ms) < e, (2-17) 

where 

.7 L 

Ve( x ) ~^2^2 c e( e , z i)M x , z j), ( 2 - 18 ) 
j=i e=o 

and 

i)i(x,z)=Y l {a')ht{k\x-z\), a' = -?—^-, z e D, xeM. 3 \D. 

\x-z\ 

(2.19) 

(ii) . 

IK - v\\ L 2 {s) < e (2.20) 

and 

\\\v e - v\\\ = 0{e) , e^0, (2.21) 



G 



where 

III ' HI = || ' \\H™ C (D') + || ' ||z,2(D';(l + |x|)-V) ) 

7 > 1 , m > is an arbitrary integer, and H m is the Sobolev space. 

Proof, (i) Note that in Theorem 2.1 we had L = L(e), while now we have L 
fixed and J = J(e). But the proof of Theorem 2.3 is similar to that of Theorem 
2.1. Let {zj}'jL 1 be a countable dense subset of U. To establish (2.17) it is 
sufficient to show that 

H :=span{M s , z j) ■ < £ < L, j = 1, 2, ...} = L 2 (S) . (2.22) 

Suppose that there exists p £ L 2 (S), p ^ such that p _L i? in L 2 (S). Define 
the single-layer potential by 

I s -3/1 

Then 

W(^) = / Vo(s, Zj)p(s) ds = (2.24) 

for j = 1,2,.... 

The continuity of the single-layer potential in M 3 implies that W(y) = for 
all y e J7. The rest of the proof is the same as in Theorem 2.1. 

□ 

Remark. Functions {^zj^Lo are linearly independent on S. Indeed, if some 
finite combination of these functions vanishes on S, then it also vanishes in the 
exterior domain D', since such a combination is an outgoing solution of the ex- 
terior Dirichlet problem with zero boundary conditions on S. In particular, such 
a combination also vanishes on Sr. Since the spherical functions are orthogonal 
on Sr, it implies that such a combination must be trivial. 

See Sections 6 and 7 for an extension of the MRC method to static problems, 
and to scattering by periodic structures, respectively. 

3 Iterative MRC algorithms 

Let z be a point in the interior of the obstacle D, and i£l 3 \D. Recall that 
Mx,z)=Y e (a')h e (k\x-z\), (3.1) 

ikr 

where hi(r) are the spherical Hankcl functions, normalized so that hi(r) ~ ^— - 
as r — > +oo. 

Noniterative MRC. 

In this MRC implementation one chooses a set of interior points H = {xj e 
D, j = 1, 2, J, J > 0} and minimizes 
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J L 

*(c) - \\u Q {s) + J2J2c eij tp e (s,Xj)\\ms), (3.2) 

j=l £=0 

over c G C^, where c = {cy}. That is, the total field w(s) = u (s) + v(s) 
is desired to be as close to zero as possible at the boundary S, to satisfy the 
required condition for soft scattering. If the resulting residual r' mm = min $ is 
smaller than the prescribed tolerance e, then the procedure is finished, and the 
sought scattered field is 

J L 

If the residual r' mm > e then the method fails. This approach, which can be 
called a Multi-point MRC, is justified by Theorem 2.3. See [13, 32, 10] for 
details and results of numerical experiments. The results show that the method 
is very efficient for domains D of a nearly spherical shape, i.e. without elongated 
parts. Clearly, the only limitation in this method is the computer resources. The 
method becomes impractical for large sets of interior points H . 

To remedy this situation one can use iterative MRC implementations, of 
which we describe the one based on a random choice of interior points, and 
another one based on an optimal choice of such points. 

Iterative MRC with a random choice of points. 

Informally, the Random Multi-point MRC algorithm can be described as 
follows. 

Fix a J > 0. Let xj,j = 1,2,..., J be a batch of points randomly chosen 
inside the obstacle D. 

Let g(s) = Uo(s), s G S, and minimize the discrepancy 

J L 

*(c) = \\g(s) + J2J2 c ^Ms,Xj)\\ms) (3-3) 
j=i e=o 

over c G C^, where c = {cej}- That is, the total field u(s) — g(s) + v(s) 
is desired to be as close to zero as possible at the boundary S, to satisfy the 
required condition for soft scattering. If the resulting residual r mm = min $ is 
smaller than the prescribed tolerance e, then the procedure is finished, and the 
sought scattered field is 

J L 

If, on the other hand, the residual r mm > e, then we continue by trying to 
improve on the already obtained fit in (3.3). Adjust the field on the boundary by 
letting g(s) := g(s) + v e (s), s G S. Create another batch of J points randomly 
chosen in the interior of D, and minimize (3.3) with this new g(s). Continue 
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with the iterations until the required tolerance e on the boundary S is attained. 
In each iteration accumulate new interior points Xj and the corresponding best 
fit coefficients C£j. After the desired tolerance is reached, the sought scattered 
field v e is computed anywhere in D'. 

Here is a precise description of the algorithm. 

Random Multi-point MRC. 

For Xj <G D, and I > functions ipe(x,Xj) are defined as in (3.1). 

1. Initialization. Fix e > 0, L > 0, J > 0, N max > 0. Let n = 0, and 

.g(s) = u ( s ), seS. 

2. Iteration. 

(a) Let n := n + 1. Randomly choose J points xj™' € D, j = 1, 2, . . . , J. 

(b) Minimize 

*(c) = \\g( 8 ) + £ ^ 4 n) )|| L2(s) 
over c G C^, where c = {c£j}. 

Let the minimum of $ be attained at c( n ) = {c^)}, j = 1,2, . . . , J, 
and the minimal value of $ be r mm . 

3. Stopping criterion. 

(a) If r mm < e, then stop. Compute the approximate scattered field 
anywhere in D' by 

n J L 

^W : =EEE c S^(^4 fe) )' ( 3 - 4 ) 

k=i j=i e=o 

(b) Ifr mm >e, and n ? N max , let 

,7 L 

:= ffW + EE ^^ 8 ' 1 ?')' x E S 
j=i e=o 

and repeat the iterative step (2). 

(c) If r mm > e, and n = N max , then the procedure failed. 

Numerical experiments based on this method are presented in the next sec- 
tion. The method is relatively slow, and it can be improved by choosing the 
interior points in some optimal way. 

Iterative MRC with an optimal choice of points. 

In this case the interior points z\, Z2, ■■■ in D are chosen one at a time, and 
their placement is not random. Rather, the discrepancy ^ is minimized not 



9 



only with respect to the coefficients c, but also with respect to the position of 
these points Zj. 

Let gi(s) = Uo(s) — u (s, a), s e S. 

Minimize 

L 

$(zi,c(zi)) := mm mm \\ gi (s) + V c t ip e {s, z)\\ L ^(S), (3.5) 
zeD c ec N j—^ 

where c = {c e } = {ce m }o<e<L,-e< m <e, L > is a fixed integer, and J2e=o '■= 

E/=oEm=-/- Lct 

L 

vi(x) = y]ce(zi)ipe(x,zi), ce(zi) = c e (zi, a). (3.6) 

The requirement (3.5) means that the total field u(s) — gi(s) + v\(s) has to 
be as close to zero as possible on the boundary S, so that it approximates best 
the Dirichlet boundary condition in (1.1). This is achieved by varying the in- 
terior point z G D and choosing the coefficients c(z) e C N giving g\ + v\ the 
best fit to zero on the boundary S. Let the minimum in (3.5) be attained at 
zi e D. If the resulting value of the residual r min = $(zi, c(zi)) is smaller than 
the prescribed tolerance e, than the procedure is finished. The sought approx- 
imate scattered field is vi(x), ief (see Theorem 2.3), and the approximate 
scattering amplitude is 

L 

A 1 (a',a)=e- ika '- Zl J2 c ^) Y ^')- (3-7) 

1=0 

Note that cn(z\) = ce(zi, a). 

The expression for Ai(a', a) in (3.7) is obtained from (3.6) by letting \x\ — ► 
cxd in x = a'\x\, because of our normalization 

ik \x | ( / ^ \ \ 

^(fcl^i—jl + O^-Jj, \x\ -» oo, (3.8) 

and \x — z\ = \x\ — a' ■ z + 0(l/\x\) as \x\ — > oo. 

If, on the other hand, the residual r mm > e, then we continue by trying to 
improve on the already obtained fit in (3.5) as follows. Adjust the field on the 
boundary by letting 52(5) = <?i(s) +i>i(s), s G S, and do the minimization (3.3) 
with 52(5) instead of gi(s), etc. Continue with the iterations until the required 
tolerance e on the boundary S is attained. At the same time keep track of 
the changing approximate scattered field v n (x), and the scattering amplitude 
A n (a',a). In this construction g n +i = u + v n on S. The goal of (3.3) is to 
obtain g n — > in L 2 (S) as n — > 00, yielding uq + w„ — > in L 2 (S). According 
to Theorem 2.3, this gives an approximate scattered solution u„ on D' to (1.1)- 
(1.2). 

Here is a precise description of the algorithm. 
MRC method with optimal choice of sources. 
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1. Initialization. Fix e > 0, L > 0, N max > 0. Let n = 0, Uo(a;) 
0, Ao(a' , a) = 0, and gi(s) = wo(s), s£ S. 

2. Iteration. 

(a) Increase the value of n by 1. 

(b) Minimize 

L 

$(z n ,c(z n )) := mm min \\g n (s) + V] c^^s, z)\\ L 2, s) , 
z£D c eC N ^ 

with the minimal value attained at z n € D, c(z n ) e C^. 

(c) Let 



v„{x) = v n --i(x) + ^2ce(z n )il)e(x,z n ), x e £>', 
i=o 

L 

A n (a', a) = A n ^(a', a) + e ^ ka '^ ^ c e (z n )Y e (a'), 



1=0 



and 



g n +i(s) = g n (s) + ^2 ce{z n )ipe{s, z n ), s e S, 
l=o 

that is g n +i(s) = uo(s) +v n (s), s e S. 
(d) Let 

r mm :=<f>(z n ,c(z n )). 

3. Stopping criterion. 

(a) If r mm < e, then stop; v n (x) is the approximate scattered field, and 
A n (a\a) is the approximate scattering amplitude. 

(b) If r mm > e, and n < N max , then repeat the iterative step (2). 

(c) If r mm > e, and n = N max , then the procedure failed. 

4 Numerical Experiments for Random Multi- 
point MRC 

In this section we describe numerical results obtained by the Random Multi- 
point MRC method for 2D and 3D obstacles. We also compare the 2D results to 
the ones obtained by the Multiple-point MRC described above, and introduced 
in [13]. The results of [13] show a favorable comparison of the Multi-point MRC 
method with the Boundary Integral Equation Method. Further improvements 
are attained with the Random Multi-point MRC method, for which the Multi- 
point MRC is just the first iteration. 
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Note that in a 2D case instead of (3.1) one has 

1 p l (x,x j ) = Hl 1 \k\x-x j \)e iW \ 

where (x — Xj)/\x — Xj\ = e l0j . 

For a numerical implementation choose M nodes {t m } on the surface S of 
the obstacle D. After the interior points Xj, j = 1, 2, J are chosen, form N 
vectors 

aW = {^(t m ,^)}^i, 

n = 1, 2, . . . , N of length M. Note that N = (2L + I) J for a 2D case, and 
iV = (L + 1) 2 J for a 3D case. It is convenient to normalize the norm in R M by 

1 M 

H b H 2 = ^ Z>™| 2 . b= (6 1 ,6 2 ,...,6 M ). 

rn—1 

Then ||uo|| - 1. 

Now let b = {g{t m )}m=ii m the Random Multi-point MRC (see Section 3), 
and minimize 

$(c) = ||b + Ac||, (4.1) 

for c <E C*, where A is the matrix containing vectors aS n \ n = 1,2, . . . , N as 
its columns. 

The Singular Value Decomposition (SVD) method (sec e.g. [24]) is used to 
minimize (4.1). According to the SVD, the matrix A is represented as 

A = UWV H , 

where the M x N matrix U has orthonormal columns u^"\ n = 1,...,N, the 
square N x N matrix V has orthonormal columns v^ n \ n = 1, . . . , N, and the 
diagonal square N x N matrix W = (w n )^ =1 is composed of the (nonnegative) 
singular values of A. 

Let P C {1, 2, . . . , N} be defined by 

P = {n : w n > w min } 

for some positive constant w min . 
Compute the normalized residual 



l — /||b|i 2 -yi<u(™),b>i 2 . 



The minimizcr c is given by 



c=V — <nW,b>vW 



n£P 



Small singular values w n < w m i n of the matrix A are used to identify and 
delete linearly dependent or almost linearly dependent combinations of vectors 
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Let r mm be the residual, i.e. the minimal value of $(c) attained after N max 
iterations of the Random Multi-point MRC method (or when it is stopped). 
For a comparison, let r^™ be the residual obtained in [13] by the Multi-point 



We have conducted 2D numerical experiments for four obstacles: two ellipses 
of different eccentricity, a kite, and a triangle. The M=720 nodes t m were uni- 
formly distributed on the interval [0, 2w], used to parametrize the boundary S. 
Each case was tested for wave numbers k = 1.0 and k = 5.0. Each obstacle was 
subjected to incident waves corresponding to a = (1.0, 0.0) and a = (0.0, 1.0). 

The results for the Random Multi-point MRC with J = 1 are shown in Table 
1, in the last column r mm . In every experiment the target residual e = 0.0001 
was obtained in under 6000 iterations, in about 2 minutes run time on a 2.8 
MHz PC. 

In [13], we have conducted numerical experiments for the same four 2D 
obstacles by a Multi-point MRC, as described in the beginning of this section. 
The interior points Xj were chosen differently in each experiment. Their choice 
is indicated in the description of each 2D experiment. The column J shows the 
number of these interior points. Values L = 5 and M = 720 were used in all 
the experiments. These results are shown in Table 1, column rJJ™. 

Thus, the Random Multi-point MRC method achieved a significant improve- 
ment over the Multi-point MRC. 

Experiment 2D-I. The boundary S is an ellipse described by 



The Multi-point MRC used J = 4 interior points xj = 0.7r( 7r(i 2 j = 
1, . . . , 4. The run time was 2 seconds. 

Experiment 2D-II. The kite-shaped boundary S (see [9], Section 3.5) is 
described by 

r(t) = (-0.65 + cost + 0.65 cos 2t, 1.5 sint), 0<t<2n. (4.3) 

The Multi-point MRC used J = 16 interior points Xj = 0.9r( 7r(j ~ 1) ), j = 
1, . . . , 16. The run time was 33 seconds. 

Experiment 2D-III. The boundary S is the triangle with vertices at 
(-1.0,0.0) and (1.0, ±1.0). The Multi-point MRC used the interior points 
Xj = 0.9r( 7r ^ g " 1 ^ ), j = 1, . . . , 16. The run time was about 30 seconds. 

Experiment 2D-IV. The boundary S is an ellipse described by 



The Multi-point MRC used J = 32 interior points xj = 0.95r( 7r<3 16 1] ), j = 
1, . . . , 32. The run time was about 140 seconds. 

The 3D numerical experiments were conducted for 3 obstacles: a sphere, 
a cube, and an ellipsoid. We used the Random Multi-point MRC with L = 



MRC. 



r(t) = (2.0 cost, sin t), < t < 2tt . 



(4.2) 



r(t) = (0.1 cost, sin t), < t < 2ir . 



(4.4) 
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Table 1: Normalized residuals attained in the numerical experiments for 2D 
obstacles by Random Multi-point MRC, ||uo|| = 1. 



Experiment 


J 


k 


a 




r mm 


I 


4 


1.0 


(1.0,0.0) 


0.000201 


0.0001 


Ellipse 


4 


1.0 


(0.0,1.0) 


0.000357 


0.0001 




4 


5.0 


(1.0,0.0) 


0.001309 


0.0001 




4 


5.0 


(0.0,1.0) 


0.007228 


0.0001 


II 


16 


1.0 


(1.0,0.0) 


0.003555 


0.0001 


Kite 


16 


1.0 


(0.0,1.0) 


0.002169 


0.0001 




16 


5.0 


(1.0,0.0) 


0.009673 


0.0001 




16 


5.0 


(0.0,1.0) 


0.007291 


0.0001 


III 


16 


1.0 


(1.0,0.0) 


0.008281 


0.0001 


Triangle 


16 


1.0 


(0.0,1.0) 


0.007523 


0.0001 




16 


5.0 


(1.0,0.0) 


0.021571 


0.0001 




16 


5.0 


(0.0,1.0) 


0.024360 


0.0001 


IV 


32 


1.0 


(1.0,0.0) 


0.006610 


0.0001 


Ellipse 


32 


1.0 


(0.0,1.0) 


0.006785 


0.0001 




32 


5.0 


(1.0,0.0) 


0.034027 


0.0001 




32 


5.0 


(0.0,1.0) 


0.040129 


0.0001 



0, Wmin = 10~ 12 , and J = 80. The number M of the points on the boundary 
S is indicated in the description of the obstacles. The scattered field for each 
obstacle was computed for two incoming directions on = (9, <p), i = 1,2, where 
<p was the polar angle. The first unit vector ot\ is denoted by (1) in Table 
2, ai = (0.0, 7r/2). The second one is denoted by (2), a 2 = (7r/2,7r/4). A 
typical number of iterations N iter and the run time on a 2.8 MHz PC are also 
shown in Table 2. For example, in experiment I with k — 5.0 it took about 
700 iterations of the Random Multi-point MRC method to achieve the target 
residual r mm — 0.001 in 7 minutes. 

Experiment 3D-I. The boundary S is the sphere of radius 1, with M = 
450. 

Experiment 3D-II. The boundary S is the surface of the cube [— 1, l] 3 
with M = 1350. 

Experiment 3D-III. The boundary S is the surface of the ellipsoid x 2 /16 + 
y 2 + z 2 = 1 with M = 450. 

In the last experiment the run time could be reduced by taking a smaller 
value for J. For example, the choice of J = 8 reduced the running time to about 
6-10 minutes. 

Numerical experiments show that the minimization results depend on the 
choice of such parameters as J, w m i ni and L. 
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Table 2: Normalized residuals attained in the numerical experiments for 3D 
obstacles by Random Multi-point MRC, ||uo|| = 1. 



Experiment 


k 




r mm 


N lter 


run time 


I 


1.0 




0.0002 


1 


1 sec 


Sphere 


5.0 




0.001 


700 


7 min 


II 


1.0 


(1) 


0.001 


800 


16 min 


Cube 


1.0 


(2) 


0.001 


200 


4 min 




5.0 


(1) 


0.0035 


2000 


40 min 




5.0 


(2) 


0.002 


2000 


40 min 


III 


1.0 


(1) 


0.001 


3600 


37 min 


Ellipsoid 


1.0 


(2) 


0.001 


3000 


31 min 




5.0 


(1) 


0.0026 


5000 


53 min 




5.0 


(2) 


0.001 


5000 


53 min 



5 Numerical Experiments for Optimal Choice 
MRC 

In this section we describe numerical results obtained by the MRC method with 
the optimal choice of sources for 2D and 3D obstacles. The notations are kept 
the same as in the previous section for the Random Multi-point MRC. As there, 
one has to minimize 

$(*,c) = ||b + Ac||, (5.1) 

in every iterative step, but, in addition, the residual is minimized with respect 
to the interior point z e D. 

There is a variety of methods to minimize &(z,c(z)), since after the mini- 
mization in the coefficients c(z) by the SVD it is just a 2D or 3D minimization 
in the region D. Our choice was the Powell's method which imitates the con- 
jugate gradients approach, but docs not require analytical expressions for the 
gradient. The Brent method was used for a line minimization, see [24, 16] for 
details. The Powell's algorithm is also described below in Section 10. 

In addition to the four obstacles considered for the Random Multi-point 
MRC, the circle \r\ = 1 was tested to check if the Optimal point MRC was 
able to find the scattered field just after one iteration, since, in this case, the 
optimal point was in the origin. The result is in Table 3, experiment number V. 
The column N iter shows the number of iterations (number of source points) at 
the end of the iterative process. The process was stopped after the algorithm 
reached the sought tolerance e = 0.002, or N max = 100. Values L = 5 and 
M = 720 were used in all 2D experiments. 

The last column (MRC - BIEM)/ BIEM shows the discrepancy in the 
scattering amplitude computed by the MRC and BIEM methods. The values 
shown are the Li norms of the difference of the scattering amplitude obtained 
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Table 3: Normalized residuals attained in the numerical experiments for 2D 
obstacles by the Optimal Choice MRC, ||uo|| = 1. 



Experiment 


k 


a 


Nuer 


r mm 


MRC-BIBM)/BIEM 


I 


1.0 


(1.0,0.0) 


20 


0.0010 


0.0001 


Ellipse 


1.0 


(0.0,1.0) 


20 


0.0018 


0.0001 




5.0 


(1.0,0.0) 


53 


0.0010 


0.0001 




5.0 


(0.0,1.0) 


45 


0.0020 


0.0001 


II 


1.0 


(1.0,0.0) 


53 


0.0020 


0.0001 


Kite 


1.0 


(0.0,1.0) 


32 


0.0020 


0.0001 




5.0 


(1.0,0.0) 


75 


0.0020 


0.0003 




5.0 


(0.0,1.0) 


68 


0.0020 


0.0001 


III 


1.0 


(1.0,0.0) 


55 


0.0020 




Triangle 


1.0 


(0.0,1.0) 


48 


0.0017 






5.0 


(1.0,0.0) 


72 


0.0019 






5.0 


(0.0,1.0) 


80 


0.0020 




IV 


1.0 


(1.0,0.0) 


100 


0.0041 


0.0008 


Ellipse 


1.0 


(0.0,1.0) 


100 


0.0027 


0.0000 




5.0 


(1.0,0.0) 


100 


0.0058 


0.0004 




5.0 


(0.0,1.0) 


100 


0.0037 


0.0012 


V 


1.0 


(1.0,0.0) 


1 


0.0000 


0.0001 


Circle 


5.0 


(1.0,0.0) 


21 


0.0020 


0.0001 



by MRC and BIEM, over the L 2 norm of the scattering amplitude obtained 
by BIEM. We followed [9] for the BIEM implementation using 64 points on 
the boundary S in every 2D experiment. No comparison is provided for a 
triangular obstacle, since it requires a complete rewriting of the BIEM code to 
accommodate the corner points. No such rewriting is required for the MRC 
method. Table III shows that for the value of tolerance e = 0.002 the computed 
scattering amplitude is in an excellent agreement with the scattering amplitude 
computed using BIEM. Results for 3D obstacles are provided in Table 4. 

Concerning the efficiency of the methods: for simple geometries the Multi- 
point MRC (see [13]) is the fastest, provided that the required accuracy can be 
achieved by a relatively small number J of the interior points (sources) used 
simultaneously. This assures the resulting matrices being of a manageable size. 
Otherwise, one has to use Random, or Optimal choice MRC, which take a 
significantly longer time to run, but can accomplish the solution of scattering 
problems untractable by single step methods, such as the Multi-point MRC 
or BIEM. While the precision of the Random-point MRC was higher in the 
conducted experiments, the optimally placed MRC method achieves an order of 
magnitude improvement in run time over it. 
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Table 4: Normalized residuals attained in the numerical experiments for 3D 
obstacles by the Optimal Choice MRC, ||uo|| = 1. 



Experiment 


k 


Qj 


Niter 


r mm 


I 


1.0 




1 


0.0000 


Sphere 


5.0 




43 


0.0019 


II 


1.0 


(1) 


12 


0.0019 


Cube 


1.0 


(2) 


7 


0.0019 




5.0 


(1) 


70 


0.0019 




5.0 


(2) 


35 


0.0020 


III 


1.0 


(1) 


12 


0.0016 


Ellipsoid 


1.0 


(2) 


35 


0.0020 




5.0 


(1) 


55 


0.0020 




5.0 


(2) 


67 


0.0020 



6 MRC for static problems 

In this Section we follow [35] and [34], Chapter 12. Consider a bounded domain 
D C M 3 with a Lipschitz boundary S, D C Br := {x : \x\ < R}. Denote the 
exterior domain by D' = M. 3 \D. Let S 2 denote the unit sphere in K 3 . Consider 
the problem: 

V 2 v = in D', v = f on S, (6.1) 

v := O ^ , r := \x\ -> oo. (6.2) 

Denote by Ye (a), a E S 2 the orthonormal spherical harmonics, Ye — Ye m , 
— £ < m < £, and let harmonic functions He(x) be defined by 

H t {x):=™, £>0, a:= X - & S 2 . 
In the region r > R the solution to (6.1)-(6.2) is: 

oo 

v(x) =^ceH t (x), r>R. (6.3) 
e=o 

The summation in (6.3) and below includes summation with respect to m, — £ < 
m < £, and ce — Q iTO are some coefficients determined by /. 

The series (6.3) in general does not converge up to the boundary S. Our 
aim is to give a formulation of an analog of the Modified Rayleigh Conjecture 
(MRC) from [27], which can be used in numerical solution of the boundary- 
value problems. The authors hope that the MRC method for static problems 
can be used as a basis for an efficient numerical algorithm for solving boundary- 
value problems for Laplace equations in domains with complicated boundaries. 
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In above sections such algorithms were developed on the basis of MRC for 
solving boundary-value problems for the Helmholtz equation. Although the 
boundary integral equation methods and finite elements methods are widely 
and successfully used for solving these problems, the method, based on MRC, 
proved to be competitive and often superior to the currently used methods. 

We discuss the Dirichlet condition but a similar argument is applicable to 
the Neumann and Robin boundary conditions. Boundary-value problems and 
scattering problems in rough domains were studied in [28] and in [34], Chapter 
9. 

Let us present the basic results on which the MRC method is based. 
Fix e > 0, an arbitrary small number. 

Lemma 6.1. There exist L = L(e) and eg — such that 

L(e) 

\\J2 c e(e)H e -f\\ L 2 {s) <e. (6.4) 
If (6.4) and the boundary condition (6.1) hold, then 

L(e) 

\\ve~ v\\ L 2 (s) < e, v e := ^2c e (e)H e . (6.5) 

Lemma 6.2. If (6.4) holds then 

\\v e - v\\ = 0(e) e^O, (6.6) 

where || • || := || • \\ H ™ C (D') + II ■ I \l 2 (d' ; (i+\x\)-'), 1 > 1, m > is an arbitrary 
integer, and H m is the Sobolev space. 
In particular, (6.6) implies 

\\ve ~ v\\ L 2 {Sr) =0(e) e^O. (6.7) 

Let us formulate an analog of the Modified Rayleigh Conjecture (MRC): 
Theorem 6.1 (MRC): For an arbitrary small e > there exist L(e) and 
ce(e),0 < £ < L(e), such that (6.4) and (6.4) hold. 
Theorem 6.1 follows from Lemmas 6.1 and 6.2. 

For the Neumann boundary condition one minimizes 1 1 d ^ ,e =° [ Ce ^ — /||x,2(s) 
with respect to q. Analogs of Lemmas 6.1-6.2 are valid and their proofs are 
essentially the same. 

If the boundary data / € C(S), then one can use C(S)— norm in (6.4)-(6.7), 
and an analog of Theorem 6.1 then follows immediately from the maximum 
principle. 

To solve problem (6.1)-(6.2) using MRC, fix a small e > and find L(e) 
and ce(e) such that (6.4) holds. This is possible by Lemma 6.1 and can be 

done numerically by minimizing || ctHt — /||l 2 (s) : ~ <?K c ij i c l)- If the 

minimum of (f> is larger than e, then increase L and repeat the minimization. 
Lemma 6.1 guarantees the existence of such L and a that the minimum is 
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less than e. Choose the smallest L for which this happens and define v e := 
^2e =0 cgHi. Then, by Lemma 6.2, v e is the approximate solution to problem 
(6.1)-(6.2) with the accuracy 0(e) in the norm || • ||. 
Proof of Lemma 6.1. We start with the claim: 

Claim: The restrictions of harmonic functions Hg on S form a total set in 
L 2 (S). 

Lemma 6.1 follows from this claim. Let us prove the claim. Assume the 
contrary. Then there is a function g ^ such that J s g(s)hg(s)ds = W > 0. 
This implies V(x) := J s g(s)\x - s^ds = OVi £ D' . Thus V = on S, 
and since AV = in D, one concludes that V = in D. Thus g = by the 
jump formula for the normal derivatives of the simple layer potential V. This 
contradiction proves the claim. Lemma 6.1 is proved. □ 

Proof of Lemma 6.2. By Green's formula one has 

w e (x)= w e (s)G N (x,s)ds, ||w £ || L 2 (s) < e, w e := v e - v. (6.8) 
Js 

Here N is the unit normal to S, pointing into D' , and G is the Dirichlet Green's 
function of the Laplacian in D': 

V 2 G = -6(x - y) in £>', G = on S, (6.9) 

G = o[^j, r^oo. (6.10) 

From (6.8) one gets (6.7) and (6.6) with Hf£ c (D')— norm immediately by the 
Cauchy inequality. Estimate (6.6) in the region B' R := R 3 \ Br follows from the 
estimate 

\G N (x,s)\ <T^< \ X \^ R - ( 6 - n ) 

In the region Br\D estimate (6.6) follows from local elliptic estimates for w f :— 
v e — v, which imply that 

IKIU 2 (b k \.d) < ce. (6.12) 

Let us recall the elliptic estimate we have used. Let D' R := Bn\D and Sr be the 
boundary of Br. Recall the elliptic estimate for the solution to homogeneous 
Laplace equation in D' R ( see [20], p. 189): 

IKIIff°-5(Dy < c[||w e || L 2 (SH) + |K|| L 2 (S) ]. (6.13) 

The estimates ||w e ||L2(s R ) = 0(e), \\w e \\ L 2^ = 0(e), and (6.13) yield (6.6). 
Lemma 6.2 is proved. □ 



7 MRC for scattering by periodic structures 

Determination of fields scattered by periodic structures is of a great impor- 
tance in modern diffractive optics, and there is a vast literature on both the 
direct and inverse problems of this type, see, for example [23]. Still, an efficient 
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computation of such fields presents certain difficulties. In the next Sections we 
present some theoretical background, a modification of the MRC method, and 
numerical results for such a scattering, see [29] . 

For simplicity we consider a 2-D setting, but our arguments can be as easily 
applied to n-dimensional problems, n > 2. Let / : R — ► R, f(x + L) = f(x) be 
an L-periodic Lipschitz continuous function, and let D be the domain 

D = {(x,y) : y > f(x), x E R}. 

Without loss of generality we assume that / > 0. If it is not, one can choose 
the origin so that this assumption is satisfied, because M := sup 0<x<i \f{x)\ < 
oo. 

Let x = (x, y) and u(x) be the total field satisfying 

(A + /c 2 )m = 0, xeD, k = const > (7.1) 

u = on S: = dD, (7.2) 

u = u +v, u : = e ika ' x , (7.3) 

where the unit vector a — (cos#, — sin0), < 9 < ir/2, and v(x) is the scattered 
field, whose asymptotic behavior as y — > oo will be specified below, and 

u(x + L,y) = vu(x,y), u x (x + L, y) = vu x {x, y) in D, v := e lkLcosd . (7 A) 

Conditions (7.4) arc the qp (quasiperiodicity) conditions. To find the 

proper radiation condition for the scattered field u(x) consider the spectral 
problem 

ip"+£ 2 ip = 0, 0<x<L, (7.5) 

ip(L) = MO), V'{L) = i/^(0) (7.6) 

arising from the separation of variables in (7.1)- (7.4). This problem has a dis- 
crete spectrum, and its eigenfunctions form a basis in L 2 (0,L). One can show 
that the corresponding eigenfunctions are e li j x and e~ lt j x with 



lj = kcos9 + ^-, or q = -k cos 9 + j = 0, ±1, ±2, . . . 

L J L 



We will use the system e 1 ^ x , which forms an orthogonal basis in L 2 (0, L). One 
has: 

f e iltx e -i^* dx = f e^O-™) dx = 0, 3+ m. 
Jo Jo 

The normalized eigenfunctions are 



e j 



<Pi{x) = -j=r, J =0,±1,±2,... 

These functions form an orthonormal basis of L 2 (0, L). 
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Let us look for v(x) = v(x, y) of the form 

oo 

v(x, y) = ^ CjVj(y)<fij(x), y > M, Cj = const. (7.7) 

j=—oo 

For y > M, equation (7.1) implies 

v'; + (e-if) V] = Q. (7.8) 

Let us assume that I 2 ^ k 2 for all j. Then 

Vj(y) = e^y, (7.9) 
where, for finitely many j, the set of which is denoted by J, one has: 

uj = (k 2 - i 2 ) 1 / 2 > 0, if I 2 < k 2 , j e J, (7.10) 

and 

Mi = ~~ k 2 ) 1 / 2 , if l 2 >k 2 ,j$J. (7.11) 

The radiation condition at infinity requires that the scattered field v{x, y) 
be representable in the form (7.7) with Vj(y) defined by (7.9)-(7.11). 

The Periodic Scattering Problem consists of finding the solution to (7.1)- 
(7.4) satisfying the radiation condition (7.7), (7.9)-(7.11). 

The existence and uniqueness for such a scattering problem is established in 
[29] . In [1] the scattering by a periodic structure was considered earlier, and was 
based on a uniqueness theorem from [12]. There are many papers on scattering 
by periodic structures, of which we mention a few [1, 2, 4, 5, 8, 17, 18], 
[22, 21, 23, 36]. The Raylcigh conjecture is discussed in several of the above pa- 
pers. It was shown (see e.g. [23, 3]) that this conjecture is incorrect, in general. 
As we have already discussed in the previous sections, the Modified Raylcigh 
Conjecture is a theorem proved in [27] for scattering by bounded obstacles. 

The main ingredient in the solution is an analog to the half-space Dirichlet 
Green's function. The function g = g(x,£,k) can be constructed analytically 
(x= (x 1 ,x 2 ),£= (£1,62)): 

ff(x,0 = ^2<Pj(xi)ipj(€i)gj(x 2 ,&,k), (7.12) 



9 3 : = 9j(x2,b,k) 



Vj(x 2 )ipj{b), x 2 >& 

Vj(&)lpj(X2), X 2 <& 



1>i = (Mj-rV^sin^-fe + b)}, a 3 = [k 2 - X 2 ] 1 ' 2 , Vj (x 2 ) = e*>*\ 
where 

+ {k 2 - fyfy = 0, ^(-b) = 0, W[ Vj ,i>j] = 1, Xj = kcos(9) + 
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and W[v, ?/>] is the Wronskian. 

The function g is analytic with respect to k on the complex plain with cuts 
along the rays Xj — it, < r < oo, j = 0, ±1, ±2, .., in particular, in the region 
3fc > 0, up to the real positive half-axis except for the set {Aj}j=o,±i,±2,...- 

Choose b > such that k 2 > is not an eigenvalue of the problem: 



(A + k 2 )ip = 0, in L>_ fc := {(x,y) : -b<y<f(x), < x < L}. (7.13) 

(7.14) 



ip\y=-b = 0, Vjv = on S, 

■0(x + L, y) = vij){x, y), ip x (x + L,y) = vip x {x, y). 



One has 



(A + /c 2 ).g = -(5(x-0, x= (x u x 2 ), ^ = (6,6), 
xe{(x,y) : — 6 < y < oo, < x < L}, 



(7.15) 



ff| w =-6 = 0. (7.16) 

Raylcigh conjectured [36] ("Rayleigh hypothesis") that the series (7.7) con- 
verges up to the boundary Sl- This conjecture is wrong ([23]) for some f(x). 
Since the Rayleigh hypothesis has been widely used for numerical solution of 
the scattering problem by physicists and engineers, and because these practi- 
tioners reported high instability of the numerical solution, and there are no 
error estimates, we propose a modification of the Rayleigh conjecture, which 
is a Theorem. This MRC (Modified Rayleigh Conjecture) can be used for a 
numerical solution of the scattering problem, and it gives an error estimate for 
this solution. Our arguments are very similar to the ones in [27]. 

Rewrite the scattering problem (7.1)- (7.4) as 

(A + k 2 )v = in D, v= -u on S L , (7.17) 

where v satisfies (7.4), and v has representation (7.7), that is, v is "outgoing", 
it satisfies the radiation condition. Fix an arbitrarily small e > 0, and assume 
that 

||u + Yl Cj(.^M<Pj(x)\\<e,0<x<L,y = f(x), (7.18) 
IjI<j(<0 

where || • || = || • \\ms L )- 

Lemma 7.1. For any e > 0, however small, and for any uo £ L 2 (Sl), there 
exists j(e) and Cj(e) such that (7.18) holds. 

Proof. Let us prove the completeness of the system {<^ J (a;)wj(/(a;))}j = o i ±i i ±2,... 
in L 2 (Sl). Assume that there is an h £ L 2 (Sl), h^0 such that 

/ h^{x)v j (f(x))ds = (7.19) 
Js L 
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for all j. From (7.19) one derives (cf. [28], p.162-163) 

V(x):=/ Mx,^ = 0,xei)_ t . (7.20) 
Js L 

Thus V = in D^, and h = ip^ ~ = 0- Lemma 7.1 is proved. □ 
Lemma 7.2. If (7.18) holds, then 

\\\v(x) - ^2 c ]( e ) v 3(y)<P]( x )\\\ < ce > Vx,y e D L , < x < L, y = /(a;), 
Ul<j(0 

w/iere c = const > does not depend on e, x, y; R > M is an arbitrary fixed 
number, and \\\w\\\ = sup x£ZAr , iR \w(x)\ + \\w\\ h1/2{Dlr) . 

See [29] for the proof. 

From Lemma 7.2 the basic result, Theorem 7.3, follows immediately: 

Theorem 7.3. MRC-Modified Rayleigh Conjecture. Fix e > 0, however 
small, and choose a positive integer p. Find 

min||u + c m( x ) v j(y)\\ : = m(p). (7.21) 
\j\<P 

Let {cj(p)} be the minimizer of (7.21). If m(p) < e, then 

«(?)= E c i(P)<PiWvM ( 7 - 22 ) 
\j\<p 

satisfies the inequality 

\\\v-v(p)\\\ <ce, (7.23) 

where c — const > does not depend on e. If m(p) > e, then there exists 
j = j(e) > p such that m(j(e)) < e. Denote Cj(j(e)) : = Cj(e) and v(j(e)) : = v e . 
Then 

\\\v - v e \\\ < ce. (7.24) 



8 Numerical solution of the periodic scattering 
problem 

According to the MRC method (Theorem 7.3), if the restriction of the incident 
field —uo(x,y) to Sl is approximated as in (7.21), then the series (7.22) ap- 
proximates the scattered field in the entire region above the profile y = f(x). 
However, a numerical method that uses (7.21) does not produce satisfactory re- 
sults as reported in [23] and elsewhere. Our own numerical experiments confirm 
this observation. A way to overcome this difficulty is to realize that the nu- 
merical approximation of the field —u \s L can be carried out by using outgoing 
solutions described below. 
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Let £ = (£1,62) € £>_&, where 6 > 0, 

!>-&:= : -&<&</(*), 0<6<L}, 

and <7(x, £) be defined as in (7.12). Then g(x, £) is an outgoing solution satisfying 
Ag + k 2 g = in Dl, according to (7.15). 

To implement the MRC method numerically one proceeds as follows: 

1. Choose the nodes Xj, i — 1,2, ...,N on the profile Sl- These points are 
used to approximate L 2 norms on Sl- 

2. Choose points C {2 \ -, £ (M) in D_ b , M < N. 

3. Form the vectors b = (w (x 4 )), and a< m ) = (g(x l7 £ (m) )), i = 1,2,...,N, 
m = 1, 2, M. Let A be the N x M matrix containing vectors sS m "> as 
its columns. 

4. Find the Singular Value Decomposition of A. Use a predetermined w m i n > 
to eliminate its small singular values. Use the decomposition to compute 

r mm =min{||b + Ac||, c e C M }, 

where 

IN| 2 = ^f>| 2 . 

i=l 

5. Stopping criterion. Let e > 0. 

(a) If r mm < e, then stop. Use the coefficients c = {ci, C2, Cm} ob- 
tained in the above minimization step to compute the scattered field 

by 

M 

v(x,y) = ^2 c m g{x,y,i (m) ). 

m—l 

(b) If r mm > e, then increase N, M by the order of 2, readjust the 
location of points G as needed, and repeat the procedure. 

We conducted numerical experiments for four different profiles. In each case 
we used L — -K,k — 1.0 and three values for the angle 9. Table 5 shows the 
resulting residuals r mm . Note that ||b|| — 1. Thus, in all the considered cases, 
the MRC method achieved 0.04% to 2% accuracy of the approximation. Other 
parameters used in the experiments were chosen as follows: TV = 256, M = 
64, w min = 10~ 8 , b = 1.2. The value of b > 0, used in the definition of 
g, was chosen experimentally, but the dependency of r mm on b was slight. 
The Singular Value Decomposition (SVD) is used in Step 4 since the vectors 
a ("») ) m — 1 ; 2,...,M may be nearly linearly dependent, which leads to an 
instability in the determination of the minimizer c. According to the SVD 
method this instability is eliminated by cutting off small singular values of the 
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Table 5: Residuals attained in the numerical experiments on MRC for periodic 
structures. 



Profile 


9 




I 


tt/4 


0.000424 




tt/3 


0.000407 




tt/2 


0.000371 


II 


tt/4 


0.001491 




tt/3 


0.001815 




tt/2 


0.002089 


III 


tt/4 


0.009623 




tt/3 


0.011903 




tt/2 


0.013828 


IV 


tt/4 


0.014398 




tt/3 


0.017648 




tt/2 


0.020451 



matrix A, see e.g. [24] for details. The cut-off value w m i n > was chosen 
experimentally. We used the truncated series (7.12) with \j\ < 120 to compute 
functions g(x,y,£). A typical run time on a 333 MHz PC was about 40s for 
each experiment. 

The following is a description of the profiles y = f(x), the nodes x, e Sl, 
and the poles € D_6 used in the computation of g(xi,^ m ^) in Step 3. For 
example, in profile I the x-coordinates of the N nodes x, <G Sl are uniformly 
distributed on the interval < x < L. The poles G £)_& were chosen as 
follows: every fourth node Xj was moved by a fixed amount —0.1 parallel to the 
y axis, so it would be within the region The location of the poles was 

chosen experimentally to give the smallest value of the residual r mm . 

Profile I. f(x) = sin(2x) for < x < L, U = iL/N, x 4 = /(**)), i = 
1,2, ...,7V, ^ = (x 4m ,y 4m -0.1), m= 1,2,..., M. 

Profile II. f(x) = sin(0.2x) for < x < L, U = iL/N, Xj = (U, f(U)), i = 
1,2,..., N, e (m) = (X4m,y4m-0.1), m= 1,2,..., M. 

Profile III. f(x) =xiov0<x< L/2, f(x) = L-x for L/2<x< L, U = 
iL/N, x, - (ti, /(*<))> * = 1, 2, N, = {x im ,y im - 0.1), m = 1, 2, M. 

Profile IV. f(x) = x for < x < L, U = 2iL/N, x f = (U,.f(ti), i = 
l,...,N/2, x, = (L,/(2(i - N/2)L/N)), i = N/2 + 1, N, ^ = (x 4m - 
0.03, j/4 m — 0.05), m = 1,2, ...,M. In this profile N/2 nodes Xj are uniformly 
distributed on its slant part, and N/2 nodes are uniformly distributed on its 
vertical portion x = L. 

The experiments show that the MRC method provides a competitive alter- 
native to other methods for the computation of fields scattered from periodic 
structures. It is fast and inexpensive. The results depend on the number of the 
internal points £( m ) and on their location. 
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Table 6: Near field values of two radiating solutions with practically the same 
far fields. 





a' 




Re v c 




Im v c 




Re v 




Im v 





.00000 


-1189 


.60834 


-227 


.35213 


-0. 


.54030 


-0 


.84147 





.31416 


-73 


.43878 


-15. 


.81270 


-0. 


.58082 


-0. 


.81403 





.62832 


1 


.94958 


0. 


.19051 


-0 


.69021 


-0 


.72361 





.94248 





.03298 


-0. 


.52343 


-0 


.83217 


-0 


.55452 


1 


.25664 


-1 


.07968 


-0. 


36021 


-0 


.95263 


-0 


.30412 


1 


.57080 


-1 


.13445 


0. 


,00027 


-1, 


.00000 


0, 


.00000 


i 


.88496 


-0 


.96294 


0. 


.31629 


-0 


.95263 


0, 


,30412 


2 


.19911 


-0 


.79021 


0. 


.55436 


-0 


.83217 


0, 


.55452 


2 


.51327 


-0 


.66472 


0. 


,71819 


-0 


.69021 


0. 


.72361 


2 


.82743 


-0 


.59154 


0, 


.81406 


-0 


.58082 





.81403 


3 


.14159 


-0 


.56768 


0. 


,84565 


-0. 


.54030 





,84147 


3 


.45575 


-0 


.59154 


0. 


.81406 


-0. 


.58082 


0, 


.81403 


3 


.76991 


-0 


.66472 


0. 


,71819 


-0. 


.69021 





.72361 


4 


.08407 


-0 


.79021 





,55436 


-0 


.83217 


0, 


.55452 


4 


.39823 


-0 


.96294 


0, 


,31629 


-0, 


.95263 


0, 


.30412 


4 


.71239 


-1 


.13445 


0. 


,00027 


-1, 


.00000 


0. 


,00000 


5 


.02655 


-1. 


.07968 


-0. 


.36021 


-0, 


.95263 


-0 


.30412 


5 


.34071 





.03298 


-0. 


.52343 


-0 


.83217 


-0, 


.55452 


5 


.65487 


1 


.94958 


0, 


.19051 


-0 


.69021 


-0 


.72361 


5 


.96903 


-73 


.43878 


-15. 


.81270 


-0, 


.58082 


-0, 


.81403 



9 Inverse scattering methods based on the MRC 

Suppose that an approximate location of the obstacle D is obtained a numerical 
inversion method, such as the Support Function Method (SFM), see [15, 30]. 
Then one can try to use the MRC method to improve the location of the bound- 
ary, see [27]. Such methods are under development by the authors, and they 
are going to be discussed elsewhere. Nevertheless, the MRC provides a tool for 
an easy construction of various examples illustrating the severe ill-posedness of 
the Inverse Scattering problem, which can be used for the algorithm's testing. 

Here we present one such example. Let the obstacle D be the unit circle 
{x e M 2 : |x| < 1}. If the incident field is u (x) = e tkx ' a , then the scattered 
field v(x) = —u (x) for x e S = dD, and its scattering amplitude is 

Aid a) = -\f^- e _i * e a{e ~ p) (9 1) 

where a' = x/|x| = e* e , and a = e tt3 . 

Let xi el 2 . Fix an integer L > 0, and let c e C 2i+1 . Form the radiating 
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solution 

L 

l=-L 

where (x — x\)/\x — x\\ = e ldl . Then its far field pattern is 



v c {x) = J2 c l Hl 1) (k\x-x 1 \)e M \ (9.2) 



AvM') = \[^ e-* (e- ika '^ £ q {-i) l e ae> j , (9.3) 

where a' = a;/|a;| = e* e . 
Fix an « 6 S 1 , and let 

r mm =mm{\\A Vc (a')- A(a',a)\\ : c e C 2i+1 } . (9.4) 

We conducted the minimization by the Singular Value Decomposition Method 
with the following values of the parameters: k = 1.0, L = 5, a = (1.0, 0.0), and 
X\ = (0.8,0.0). The L 2 norm in (9.4) was computed over M = 120 directions 
a' m uniformly distributed in the unit circle S 1 , and then normalized by VM, 
so that the identity function would have the norm equal to 1. The resulting 
value of the residual r mm = 0.00009776 indicates that the far field A{a! , a) was 
practically perfectly fit by the radiating solution of the form (9.2). However, as 
the Table 6 shows, the restrictions of the exact scattered field v, and the fitted 
field v c to the boundary S of the obstacle D are vastly different. The columns 
in Table 6 correspond to the real and the imaginary parts of the scattered fields, 
and the rows correspond to different values of the angle a' . Thus, one has to 
conclude that, as expected, a coincidence of the radiating solutions at the far 
field does not imply that the near fields are also coincidental, see [30]. 



10 Stability Index Method 

Various algorithms for direct and inverse scattering problems require global 
minimization of functions of many variables, see [30]. Since most objective 
functions contain many local minima, this is a highly nontrivial task. In several 
papers, starting with [14], the authors developed and tested the Stability Index 
Method (SIM) for global minimization. In our presentation here we follow [11], 
which also contains a convergence analysis and additional numerical results. 

The Stability Index Method combines stochastic and deterministic algo- 
rithms to find global minima of multidimensional functions. The functions may 
be nonsmooth and may have multiple local minima. The method examines the 
change of the diameters of the minimizing sets for its stopping criterion. At 
first, the algorithm uses the uniform random distribution in the admissible set. 
Then normal random distributions of decreasing variation are used to focus on 
probable global minimizers. To test the method, we have applied it to standard 
test functions of several variables. The computational results show that the SIM 
is efficient, reliable and robust. 
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Given a function / : A — ► K, our goal is to minimize it over an admissible set 
A assumed to be a bounded set in a metric space X. Typically, the structure 
of the objective function / is quite complicated. In particular, it can have 
many local minima and a non unique global minimum. To better understand 
the structure of the minima, let us introduce the minimizing sets S e of /. Let 
m = inf{/(x) : x G A}. Given an e > define 

S* e = {x e A : f(x) < m + e}, (10.1) 

or 

S e = {x e A : f(x) < f(x p ) + e}, (10.2) 

if the problem admits a global minimizer x p E A. 

Definition. Given an e > 0, let D e be the diameter of the minimizing set 
S € . We call D e the Stability Index of the minimization problem (10.1). 

We are interested in the behavior of D f as e — > 0. So, one can say that the 
problem (10.1) possesses a set of Stability Indices {D e : e > 0}, and the above 
definition should be understood in this sense. 

One would expect to obtain a stable identification for minimization problems 
with small (relative to the admissible set) stability indices. Minimization prob- 
lems with large stability indices either have distinct global minimizcrs, or the 
function / is nearly flat in a neighborhood of the global minimizer x p . In this 
situation, and with no additional information known, one has an uncertainty of 
the minimizer's choice. The stability index provides a quantitative measure of 
this uncertainty or instability of the minimization. 

In a practical minimization problem one constructs a sequence of minimizcrs 
{xi,X2, ■ •■} C A, and makes a decision when to terminate the iterations accord- 
ing to a stopping criterion. We assert that the knowledge of the Stability Index 
provides a valuable tool for the formulation of such a stopping criterion. 

Originally, we have applied the Stability Index minimization method to in- 
verse scattering problems arising in quantum mechanical scattering, [14]. Such 
potential scattering problems are important in quantum mechanics, where they 
appear in the context of scattering of particles bombarding an atom nucleus. 
One is interested in reconstructing the scattering potential from the results of 
a scattering experiment. Assuming a particular structure of the potential, the 
scattering results can be computed and compared to the given scattering data. 
Thus the inverse scattering problem is reduced to the minimization of the dis- 
crepancy (best fit to data), see [14, 31] for details. 

The goal of the SIM algorithm is to find a minimizing set S e that fits within 
a small portion of the computational domain A C M. N . Practically, we assume 
that A = [-M, M] N C R N , for an M > 0. If it is desirable to introduce different 
scales for the variables, then the algorithm should be modified accordingly. 

Let < 5 < 1. The minimization is stable if, given a global minimizer x p , 
we are able to find a minimizing set S e C C[x p ,8], where C[x p ,6] is the cube 
centered at x p € A with the side equal to 25 M. 

The next step is to define a sequence of normal distributions T n with the 
variances \i n — > 0, as n — > oo. Thus we fix an < a < 1, and let [i n — a™, n — 
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1,2,... 

Initially, for n = 0, let the To be the uniform random distribution in A. A 
special algorithm SMS, described below, determines a finite Stable Minimizing 
Set So C A. Let x € So be the minimizer in So, that is 

f(x ) = min{/(x) : z e S } . (10.3) 

If So C C[xo, <5], then the minimization is stable and the global minimizer x p = 
x . 

If, on the other hand, the above inclusion is not achieved, then one continues 
with another application of the SMS, this time using the normal distribution 
T\ with the mean at xq, and the variance fii, etc. The iterations continue until 
either S„ C C[x n ,S\ or 3fi n < 26 M. The last condition is needed to prevent 
all the trial points to be chosen too close to x n , thus preventing a reasonable 
estimate for the diameter of S n . 

Stability Index Method (SIM) 

Fix < a, S < 1. Suppose that A = [-M, M] N . 

1. Initialization. Let n = 0. Use the SMS algorithm with the uniform 
random distribution T in A to determine the minimizing set So C A and 
the minimizer x € So- Go to the Stopping Criterion (step 3) to check if 
additional iterations are needed. 

2. (n— th iteration). Let \i n — a n . Use the SMS algorithm with the normal 
random distribution T n with the mean at x„_i and the variance ji n to 
determine the minimizing set S„ci and the minimizer x n <E S n . 

3. Stopping criterion. Let C[x n ,S] be the cube centered at x n € A with the 
side equal to 28 M. 

If S„ C C[x„,<5], then stop. The minimization is stable. The estimated 
global minimizer x p is x n . 

If S n (£ C[x n , S] and 3fi n < 2SM, then stop. The minimization is unstable. 
The diameter (Stability Index) D n of S„ is a measure of the instability of 
the minimization. 

Otherwise, increase n by 1, and return to Step 2 to do another iteration. 

Note that the obtained point x p is an estimated global minimizer. 

The main part of the Stability Index Method is the SMS algorithm which 
determines stable minimizing sets S n , corresponding to the random distributions 
T n . These distributions are either uniform in A or normal with a given variation 

The SMS algorithm is an iterative algorithm. It can be called an Iterative 
Reduced Random Search method. Choose an integer K > assuming that K 
random points in S e are sufficient to estimate its diameter D f . If n > 1, then 
the calling algorithm SIM provides the minimizing set S„_i, its minimizer x n -\, 
and the variance /x„. 
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Let a batch H 1 C A oi L > K trial points be generated in the admissible 
set A according to the random distribution T n . If n = 0, then To is just the 
uniform random distribution in A. If n > 1, then T„ is the normal distribution 
with the variance /x n , and the mean at Let be the subset of K points 

from H 1 where the objective function / attains its K smallest values. That is 

max{7(wi) : u { e Q\j} < min{/(uj) : m g Q\j}. (10.4) 

Use each point Ui £ Q\j as the initial guess for a Local Minimization Method 
(LMM) of your choice, e.g. the conjugate gradient method, etc. The specific 
LMM used by us is described below. While the use of a local minimization is 
not, strictly speaking, necessary for the SIM, it provides a significant improve- 
ment in the performance of the algorithm, and is highly recommended. Thus for 
each starting point U{ <G Q\, the LMM produces a minimizer Vi € A. Let Q v 
be the set of all such minimizers. Let Q 1 be the subset of Q\j U Q v containing 
K points with the smallest values of /, and q 1 be the minimizer in Q 1 . Define 
the radius of Q 1 by 

= max{||z, -q 1 !! :z t eQ\ i = 1, 2, K}. (10.5) 

The idea of the Stability Index Method is to iteratively construct subsets Q 3 
until their diameters are stabilized. Practically, one can achieve the same goal 
by estimating and examining the radius of the set Q 3 . This also requires 
less computational effort. 

To construct the next set Q 2 generate another batch H 2 C A of L trial 
points according to the uniform random distribution, if n = 0, or, for n > 1, 
according to the normal distribution T n with the variance and the mean 
at q 1 . Let Qfj be the subset of K points from H 2 U Q 1 having the smallest K 
values of /. Apply the LMM to produce the set of minimizers Qy. Of course, 
if some point Ui € Q\j has already been used as an initial guess for the LMM 
in the previous iteration, it is excluded from the LMM application. Let Q 2 be 
the subset of Qfj U Qy containing K points with the smallest values of /. Let 
q 2 be the minimizer in Q 2 , and = max{||zi — q 2 \\ : Zi G Q 2 , i = 1,2, 
be its radius, etc. 

This way one produces a sequence of the minimizing sets Q J , j = 1,2, .... 
Let < 7 < 1, and P be a positive integer. The iterations are terminated if 
the maximum number of iterations N max is exceeded or the following Stopping 
Criterion is satisfied: 



p ^ 



i=j-p+i 



<2 7 M. (10.6) 



In either case, when the last iteration j is determined from (10.6) or j 
N max , we let S n = Q J and x n = q 3 . 

Stable Minimizing Set (SMS) algorithm 
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Fix < 7 < 1, and integer K,L > K, P, N max . Constant M, normal random 
distribution T n , its variance jj, n (for n > 1), the minimizing set S n -i, and the 
minimizer x n -\ are supplied by the calling algorithm SIM. 

1. Initialization. Let j = 1. 

• For n = 0. Generate a batch H 1 of L trial points in A C R N using 
the uniform random distribution. Let Q\j be the subset of K points 
from H 1 where the objective function / attains its K smallest values. 
Go to step 4. 

• For n > 1. Generate a batch H 1 of L trial points in A C R N using 
the normal distribution T n with the variance \i n and the mean at 
x n -\. Let be the subset of K points from H 1 U S n -i where the 
objective function / attains its K smallest values. Go to step 4. 

2. Iterative step {j > 2). 

• For n — 0. Generate a batch £P of L trial points in A C l w using 
the uniform random distribution. 

• For n > 1. Generate a batch of L trial points in ^4 C l w using 
the normal distribution T n with the variance /i„ and the mean at 

3. Let Qy be the subset of K points from W U Q^ 1 where the objective 
function / attains its K smallest values. 

4. Local minimization. Use each unflagged point Ui e Q\j as the initial guess 
for a Local Minimization Method (LMM). Let Vi e A be the resulting 
minimizer. Let Q 3 V be the set of all such minimizcrs resulting from the 
application of LMM to Q-jj. Flag all points in and Qy. 

5. Let Qi be the subset of Q\j U Q v containing K points with the smallest 
values of / and qi be the minimizer in Qi . Define the radius of Q J by 

= max{||zi- q j \\ : m e Q 1 , i = 1, 2, K}. 

6. Stopping criterion. 

• If j < P, increase j by 1 and return to step 2 for another iteration. 

• If j > f) compute the average radius during the last P iterations: 

i= 3 -P+l 

• Termination. If {R^ - R a \ < 2 7 M, or j > N maxi let S n = , 
x n = qi and exit the procedure. 
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• Otherwise, increase j by 1 and return to step 2 for another iteration. 

The SMS implementation involves a combination of stochastic (global) and 
deterministic (local) minimization methods. Generally, local searches offer more 
precision and speed than their global counterparts, so that adding a local step 
to a global minimization algorithm should yield improvement in both areas. 
Likewise, by itself, a local minimization method will very often produce points 
of considerable distance from the actual global minimizcr, that is it would be 
trapped in one of many local minima of the objective function /. Adding a 
global step helps the algorithm escape from local minima, and explore the entire 
admissible set A. The use of various normal distributions of decreasing variance 
is similar to ideas of the simulated annealing method [19]. 

The particular Local Minimization Method (LMM) used in the numerical 
experiments was a modification of Powell's minimization method in M. N , [7]. It 
was chosen with applications in mind, for which the objective function / does 
not have a convenient expression for its gradient. Either a Golden Search or 
Brent's method can be used for one dimensional minimizations, [24]. 

Modified Powell's Method 

1. Choose the set of directions Ui , i = 1, 2, . . . , N, to be the standard basis 
in R N 

u t = (0,0,...,1,...,0), 
where 1 is in the i-th place. 

2. Save the starting point pq. 

3. For i = 1, . . . , N move from pi-i along the direction Ui and find the point 
of minimum pi . 

4. Set v — pn — po- 

5. Move from p along the direction v and find the minimum. Call it p 
again. It replaces po from step 2. 

6. Repeat the above steps until a stopping criterion is satisfied. The resulting 

point is Prrnn- 

Note that f(p m in) < f(po) f° r an Y objective function / used in the Local Min- 
imization Method. 

11 Numerical results for SIM 

The Stability Index Method described in the previous sections was tested on 
several functions designed to test and compare various minimization algorithms, 
see [11] for additional test functions results. The experiments were conducted 
on a 2.8 GHz PC with 256 MB RAM. 
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In all the numerical experiments we used the same parameter values: a = 
0.8,(5 = 0.001,7 = 0.001, K = 30, Z, = 5000, P = 6, and N max = 30. For each 
test function the admissible set A is a cube [— M, M] N provided in the function's 
description together with its global minimizer. 

Test Function 1 



f(x, y) = C0S I(* + 1 ) x + *] ) I cos[(i + l)y + i] 



+ 0.5((x + 1.4213) 2 + (y + 0.80032) 2 ) 



The minimum is sought on [—5, 5] x [—5, 5]. This function has a global minimum 
at (-1.42513, -0.80032) with a function value of -186.73091, [37]. 

Test Function 2 



f(x,y) = e sin(50x) + sin(60e J/ ) + sin(70sina;) + sin(sin(80y)) 

-sin(10(x + y)) + (x 2 + y 2 )/4. 

The minimum is sought on [— 1, l] 2 . According to [6] the minimum occurs at 
approximately (-0.0244031,0.2106124) with a function value of -3.30686865. 

Test Function 3 



f{x) = - ^lOsin 2 ^) + ((W - !) 2 (! + 10sin 2 (7r2/ 4 + 1)) + (y N - 1) 2 J 

where x = [x x ,x 2 , ...,x N ) G R N , y t = 1 + 0.25(x 4 - 1), i = 1, 2, N. The 
minimum is sought on [—10,10]^. This function has a global minimum at 
x = (1, 1, 1) with a function value of 0, [37]. 

The results of the minimization using the SIM for these test functions are 
shown in Table 7. The algorithm was run 20 times for each function. It found 
the correct global minimum most of the time. The "success rate" column in 
Table 7 shows the percentage of trials in which the global minimum was found 
exactly. The " Function evaluation" column shows the average number of times 
the objective function was evaluated. Finally, Table 7 shows the average run 
time, in seconds, for a single trial run. 

12 Conclusions 

Let D be a 2D or 3D obstacle, S be its boundary, and zto be the incident field. 
Raylcigh conjectured that the acoustic field u in the exterior of the obstacle is 
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Table 7: Results of the computational experiments for SIM. 



Function 


Dimension 


Actual 


Found 


Success 


Average run 




N 


minimum 


minimum 


rate (%) 


time (seconds) 


1 


2 


-186.731 


-186.731 


100 


2 


2 


2 


-3.30687 


-3.30687 


100 


2 


3 


5 


0.00000 


0.00000 


100 


7 


3 


10 


0.00000 


0.00000 


100 


16 


3 


20 


0.00000 


0.00000 


100 


50 



given by 

CO 

u(x,a) =e ika ' x + J2M<x)i>t, ^e-=Ye{a')h e (kr), a' = -. (12.1) 
t=o r 

While this conjecture (RC) is false for many obstacles, it has been modified to 
obtain a representation for the solution of (1.1)-(1.2) and to obtain its error. 

It is proved that if v e is an outgoing solution of the Helmholtz equation in the 
exterior domain D' and uo + v t approximates zero in L 2 (S)— norm on the bound- 
ary S, then v e approximates the exact scattered field v in D', see Theorems 2.1 
and 2.3. The Modified Rayleigh Conjecture approach to obstacle scattering 
problems is based on the following observation: the functions tpe(x,z), z e D 
and their linear combinations are outgoing solutions to the Helmholtz equation 
in the exterior domain. Therefore, one just needs to find a combination of such 
functions that gives the best fit to — u on the boundary S. Then this com- 
bination approximates the scattered field everywhere in the exterior D' of the 
obstacle D and the error of this approximation is given in Theorems 2.1 and 
2.3. 

In this paper we describe several implementations of the MRC method which 
give an efficient approach to solving obstacle scattering problems for 2D and 
3D problems with complicated geometries. Our implementations of the MRC 
method worked more efficiently than the BIEM method. 

Various methods for solution of direct and inverse scattering problems re- 
quire a global minimization of the objective function. We developed the Stability 
Index Method which is a robust and efficient algorithm for global minimization. 
Its efficiency comes from a combined use of global and local minimization. The 
global (stochastic) part employs uniform and normal random distributions. It 
can be combined with local (deterministic) methods appropriate for the objec- 
tive function. The diameters of the minimizing sets (Stability Index) are used 
for a self-contained stopping criterion. The computational experiments show 
that the method was successful for various standard test functions over mul- 
tidimensional domains. No adjustment of parameters was needed in different 
tests. The method is well suited for low dimensional minimization problems. Its 
performance deteriorates for higher dimensional problems. The Stability Index 
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Method is a valuable addition to already existing global minimization methods. 
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